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FIGURES 


Figure  1.  A  plot  of  density  verses  temperature  of  water  is  shown . 5 

Figure  2.  The  geometry  for  the  two-dimensional  problem  is  shown.  The  laser  beam 
(with  its  e_1  radius)  is  illustrated  by  the  solid  circle.  All  boundaries  are  thermally 
insulated,  and  the  fluid  velocity  at  all  boundaries  is  zero  (the  non-slip  boundary 
condition  applies) . 10 

Figure  3.  Results  for  the  two-dimensional  problem  are  shown.  The  temperature 
profile  is  given  at  time  equals  one  second.  Minimum  temperature  equals  293  K; 
maximum  temperature  equals  294.6  K.  The  horizontal  centerline  is  drawn  to  indicate 
the  slight  vertical  asymmetry  of  the  temperature  distribution;  this  illustrates  an  early 
stage  of  convective  flow . 10 

Figure  4.  Results  for  the  two-dimensional  problem  are  shown.  The  temperature 
profile  at  time  equals  three  seconds  is  given.  Minimum  temperature  equals  293  K; 
maximum  temperature  equals  299  K.  The  horizontal  centerline  is  drawn  to  indicate  the 
vertical  asymmetry  of  the  temperature  distribution;  this  illustrates  effects  of  convective 
flow . 11 

Figure  5.  Results  for  the  two-dimensional  problem  are  given.  Convective  fluid 

velocity  vectors  at  time  equals  three  seconds  are  shown.  Note  the  convective 

flow  profile.  The  maximum  amplitude  of  the  velocity  is  4.4  mm/sec . 1 1 

Figure  6.  Results  for  the  two-dimensional  problem  are  given.  Temperature  profiles 
along  the  central  vertical  axis  (where  x  =  0.01  m)  are  shown  at  various  times.  The 
asymmetry  (about  the  center-y-position  at  y  =  0.01  m)  of  the  temperature  in  the 
y-direction  is  indicative  of  convective  flow.  This  asymmetric  pattern  is  shown  to 
increase  with  time . 12 

Figure  7.  Results  for  the  two-dimensional  problem  are  given.  Plots  are  shown  of 
temperature  profiles  along  the  central  horizontal  axis  (where  y  =  0.01  m)  at  various 
times.  These  results  are  presented  to  demonstrate  the  extent  of  the  right-left 
symmetry  (about  the  central-x  position  at  x=  0.01  m)  of  the  result . 12 

Figure  8.  The  geometry  for  the  three-dimensional  problems  1  and  2  is  illustrated.  The 
arrow  shows  the  direction  of  propagation  and  attenuation  of  the  laser  beam.  The 
coordinates  are  shown  in  meters,  where 

0<x< 0.02m,  0<y< 0.02m,  and  0<z<  0.04m . 13 

Figure  9.  Results  for  the  three-dimensional  problem  1  are  given.  The  temperature 
profile  at  the  front  face  is  shown  at  time  equals  four  seconds.  A  horizontal  centerline  is 
drawn  to  illustrative  vertical  asymmetry  caused  by  convection.  A  vertical  centerline 
is  drawn  to  evaluate  the  expected  symmetry  in  the  horizontal  direction.  Maximum 
temperature  equals  302.7  K;  minimum  temperature  equals  293.15  K . 14 
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Figure  10.  Results  for  the  three-dimensional  problem  1  are  given.  The  temperature 
profile  through  a  central  y-z  plane  (where  x=0.01m)  at  time  equals  four  seconds  is 
presented.  Note  the  vertical  asymmetry  of  the  temperature  profile.  Maximum 
temperature  equals  302.8K;  minimum  temperature  equals  293. 15K . 15 

Figure  11.  Results  for  the  three-dimensional  problem  1  are  given.  Temperature  profiles 
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central  y  position  (at  y=  0.01  m)  with  increasing  time.  A  small  dip  in  temperature  below 
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three  seconds.  Results  are  presented  for  various  values  of  z.  Note  the  asymmetry 
in  the  y-direction;  this  illustrates  effects  of  vertical  convection.  A  small  dip  in 
temperature  below  293K  is  caused  by  a  small  computational  inaccuracy . 16 
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Figure  14.  Results  for  three-dimensional  problem  2  are  given.  The  temperature 
profile  at  the  front  face  at  time  equals  four  seconds  is  shown.  A  horizontal 
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Figure  15.  Results  for  three-dimensional  problem  2  are  given.  The  temperature  profile 
through  the  central  vertical  plane  (at  x  =  0.01m)  at  time  equals  four  seconds  is 
presented.  The  horizontal  line  is  drawn  to  show  vertical  asymmetry  due  to  convection. 


Maximum  temperature  equals  296.2K;  minimum  temperature  equals  293. 15K . 20 

Figure  16.  Results  for  three-dimensional  problem  2  are  given.  The  temperature 
profiles  along  the  vertical  central  axis  (at  x  =  0.01m)  at  the  front  face  (where  z  =  0) 
at  various  times  are  presented . 21 


Figure  17.  Results  for  three-dimensional  problem  2  are  given.  The  temperature 
profiles,  as  a  function  of  y,  along  the  central  vertical  axis  (where  x  =  0.01)  at  time 
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Figure  18.  Results  for  three-dimensional  problem  2  are  given.  Temperature  profiles 
as  a  function  of  x  along  the  central  horizontal  axis  (y  =  0.01)  at  time  equals  four 
seconds  are  presented.  Results  are  shown  for  various  values  of  z.  These  results 
were  plotted  as  a  check  for  the  evaluation  of  the  horizontal  symmetry  of  the 
solution . 22 

Figure  19.  The  geometry  for  the  simulation  involving  the  evaporative  and 
convective  heat  transfer  at  the  z  =  0  surface  is  illustrated.  The  height  of  the 
cylinder  equals  1  cm,  and  the  radius  of  the  cylinder  equals  0.5  cm.  The  laser 
beam  enters  perpendicularly  at  the  top  (at  z  =  0)  surface.  The  problem  is  circularly 
symmetric . 24 

Figure  20.  The  simulated  temperature  profile  in  a  cylindrical  container  of  water  is 
shown.  Evaporative  cooling  and  convective  heart  exchange  take  place  at  the  top 
boundary,  where  z  equals  zero.  Results  for  time  equals  8  seconds  are  presented.  Laser 
input  power  equals  1W.  Maximum  temperature  equals  312K;  minimum  temperature 


equals  293. 15K . 27 

Figure  21.  Surface  tension  vs.  Temperature  of  water.  The  value  of  the  parameter  y 
can  be  found  by  taking  the  derivative  of  this  curve . 28 
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1. 


INTRODUCTION 


Mathematical  modeling  and  experimental  methods  have  been  used  by  researchers  to  understand 
the  optical  effects  of  laser  beam  penetration  in  the  aqueous  humor.  These  models  and 
experiments  reveal  the  phenomena  of  thermal  lensing  of  the  laser  beam  [Vincelette,  et.  al,  2007], 

Water  is  heated  by  the  power  absorbed  from  an  attenuated  laser  beam.  The  temperature  changes 
of  the  heated  water  cause  localized  refractive-index  changes  in  the  propagating  medium. 
Therefore,  a  quantitative  analysis  of  local  temperature  changes  caused  by  the  laser  heating  could 
facilitate  a  better  understanding  of  the  laser  lensing  effects. 

Laser  heating  is  caused  by  coupled  convective  and  conductive  transport  of  thermal  energy.  The 
experimental  results  of  Vincelette,  et.  al.  [2007]  suggests  the  importance  of  thermal  convection 
upon  the  optical  properties  of  laser  beam  penetration. 

When  conductive  effects  predominate,  and  convective  effects  are  negligible,  the  thermal  changes 
can  be  observed  as  circularly-symmetric  temperature  profiles.  Circularly  asymmetric 
temperature  patterns  can  be  observed  when  the  convective  effects  become  relevant. 

In  this  report,  simulation  results  for  laser  heating  of  water  will  be  presented,  and  the  conductive 
and  convective  heating  effects  will  be  evaluated. 
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2.  SIMULATION  OF  CONDUCTIVE  AND  CONVECTIVE  HEAT  TRANSFER  AND 
TEMPERATURE  CHANGES  DUE  TO  LASER  HEATING  OF  WATER: 
MATHEMATICAL  AND  THEORETICAL  BACKGROUND 

2.1  BASIC  PRINCIPLES  AND  MATHEMATICAL  EQUATIONS 


Thermal  heating  is  modeled  using  the  thermal  transport  equation  coupled  with  the  Navier-Stokes 
and  continuity  equations.  The  thermal  transport  equation  describes  temperature  changes  caused 
by  thermal  convection  and  conduction,  while  the  Navier-Stokes  and  continuity  equations 
describe  the  convective  velocity  caused  by  spatial  variations  of  the  buoyancy  force.  The  spatial 
variations  in  the  buoyancy  force  result  from  temperature-dependent  density  gradients.  Thus,  the 
Navier-Stokes,  continuity,  and  thermal  transport  equations  form  a  set  of  coupled,  partial  - 
differential  equations  from  which  the  fluid  convective  transport  and  the  temperature  variations 
can  be  determined. 


The  fluid  flow  is  governed  by  the  Navier-Stokes  equations  [Landau  and  Lifshitz,1959; 
Batchelor,  1967]: 

p  —  +  p(V»  V)V  +  VP  -  nV2V  =  F, 

dt  (i) 

and  the  continuity  condition  for  incompressible  flow 

V  •  V  =  0,  (2) 


where  p  denotes  the  density,  and  q  denotes  the  dynamic  viscosity  of  the  liquid,  t  denotes  time 
and  F  denotes  the  body  force  caused  by  spatially  varying  buoyancy  effects. 


The  equation  for  conductive  and  convective  heat  transfer  is 


pC 


'  dr 

v  dt 


\ 

+  V  •  VP 

J 


v  »(kVT)=  Qs, 


(3) 


where  C  is  the  specific  heat  (J  /  kg»K)  and  k  is  the  thermal  conductivity  of  the  liquid  (W/(m»K)). 
Qg  is  a  heat-source  term  (W/m3),  and  T  is  the  temperature  (K)  [Incropera,  et  al,  2007], 


The  heat-source  term,  Q$ ,  denotes  the  heat  absorbed  (per  unit  volume)  by  the  laser.  This  term 
can  be  expressed  as  [Torres  et  al,  1993]. 


Qs{r,z,t)  =  pa  S(r,z,t), 


(4) 
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where  /Ja  is  the  absorption  coefficient  of  the  material  (m'1).  Using  Beer’s  law  as  a  first 
approximation,  the  magnitude  of  the  laser  energy  flux  is  given  as 


P 


f 


7rR 


exp 


w 


0 r-r o) 


2  ^ 


R 


w 


exp (-juaz)f(t)  ■ 


(5) 


Here,  P  is  the  power  of  the  incident  laser  beam*  (W),  Rw  is  the  e'1  laser  beam  radius  (m),  r  is  the 
radial  coordinate  and  z  is  the  axial  coordinate.  The  position  of  the  central  axis  of  the  laser  beam 

is  denoted  as  To.  When  rectangular  coordinates  (x,y,z)  are  used,  with  a  center  axis  at  (xo,yo),  r  — 
fo  is  given  as 


r-r0  =  J(x-x0)2+(y-y0)2 


(5a) 


To  model  transient  effects  from  a  continuous- wave  laser,  f(t)  is  a  unit-step  function. 

The  Boussinesq  approximation,  which  is  explained  in  the  next  paragraph,  is  used  to  determine 
the  effect  of  temperature  gradients  upon  the  convective  velocity  of  the  liquid.  The  temperature- 
dependent  density  variations  produce  spatial  variations  of  the  buoyant  force  within  the  liquid. 
Thus,  the  x  and  y  components  of  the  gravitational  body-force,  F,  are 


Fx-0, 

Fy  =  -g  P(T) 


(6) 


where  g  is  the  acceleration  of  gravity  (9.81  m/s2),  p(T)  is  the  temperature-dependent  density  of 
the  liquid  (kg/m3).  A  plot  of  density,  p(T)  verses  temperature,  T.  for  water  is  shown  in  Figure  1 
[Comsol  Lib.,  2006]. 


To  determine  the  incident  laser  power  P,  the  reflection  losses  at  the  front  face  must  be  taken  into  account.  Thus,  P 
is  equal  to  the  power  of  the  laser  multiplied  by  the  fraction  of  the  laser  power  that  is  transmitted  through  the  front 
face  of  the  container. 
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Figure  1.  A  plot  of  density  versus  temperature  of  water  is  shown. 


Equation  (6)  is  substituted  into  the  body- force  term  at  the  right  of  the  equal  sign  in  equation  (1). 
Note  that,  as  indicated  by  equation  (6),  the  body- force  term  in  the  Navier-Stokes  equation  (1)  is 
temperature  dependent.  Also,  the  convective  term  in  the  heat-transport  equation 
(3)  is  velocity  dependent.  Thus,  equations  (1),  (2)  and  (3)  form  a  coupled  set  of  partial 
differential  equations. 

The  Boussinesq  approximation  is  used  to  determine  the  effect  of  temperature  gradients  upon  the 
convective  velocity  of  the  liquid  [Acheson,  1990;  Drazin  and  Reid,  1981;  Incropera,  et  al,  2007; 
Tritton,  pp.  188-195,  1998].  According  to  the  Boussinesq  approximation,  all  density  parameters 
in  the  equations  (1)  and  (3)  are  constants  (where  p  =  100  0kg/m3),  with  the  exception  of  the 
density  variations  involved  with  the  buoyant  force,  which  is  included  as  a  body-force  term  in  the 
Navier-Stokes  equation.  The  components  of  the  temperature-dependent  body-force  are  given  by 
equation  (6).  The  incompressibility  condition,  given  by  equation  (2),  is  included  as  part  of  the 
Boussinesq  approximation. 


2.2  DIMENSIONLESS  PARAMETERS: 


Dimensionless  parameters  involved  in  convective  and  conductive  heat  and  mass  transport  of  an 
incompressible  fluid  are  the  Reynolds  number  and  the  Rayleigh  number.  The  Reynolds  number, 
Re,  is  expressed  as 


Re  = 


Vpd 
l 1 


(7) 


Where  V  denotes  the  velocity  scale,  and  d  denotes  the  length  scale  (equal  to  a  length  of  a  side 
of  the  container).  Also,  p  denotes  the  density,  and  r\  denotes  the  dynamic  viscosity  of  the  liquid. 
The  Reynolds  number  expresses  the  order-of-magnitude  estimate  of  the  ratio  of  the  inertia  terms 
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to  the  viscous  term  in  the  Navier-Stokes  equation.  High  Reynolds  numbers  indicate  flow 
instability  for  transition  to  turbulence  [Landau  and  Lifshitz,1959].  In  our  investigation,  the 
calculated  Reynolds  numbers  are  on  the  order  of  1  -  10,  indicating  the  non-existence  of 
turbulence,  but  also  indicating  that  the  non-linear  term  and  the  time-dependent  term  in  the 
Navier-Stokes  equation  cannot  be  neglected. 

The  Rayleigh  number,  Ra,  given  as 


Ra  =  gaiA^T^P 

T]K 


(8) 


where  g  is  the  acceleration  of  gravity  (9.81  m/s2),  a  is  the  thermal  expansion  coefficient  of  the 
liquid  (K1),  AT  is  the  magnitude  of  the  spatial  variation  of  temperature  (K),  and  dx  denotes  the 
length  scale  over  which  the  thermal  gradients  are  effective.  The  thermal  diffusitivity  of  the  fluid, 
k  (m2/s),  is  given  as 


K  = 


k 

~pC 


(9) 


The  Rayleigh  number  indicates  the  instability  of  a  heated  liquid  to  the  formation  of  convective 
flow.  That  is,  at  high  Rayleigh  numbers,  convective  flow  is  expected  to  occur.  (For  example,  for 
flow  between  two  rigid  horizontal  plates,  separated  by  a  vertical  distance  dx,  instabilities  leading 
to  convective  flow  occur  when  Ra  exceeds  a  critical  Rayleigh  number  of  1708  [Velarde,  M.G., 
and  Normand,  1980;  Drazin  and  Reid,  1981].  A  critical  value  of  the  Rayleigh  number 
corresponding  to  the  physics  and  geometry  of  the  problems  considered  in  this  report  is  not 
known.)  As  indicated  by  equations  (8)  and  (9),  an  increase  in  the  dynamic  viscosity  and/or  in  the 
thermal  diffusitivity  of  the  fluid  would  inhibit  the  formation  of  convective  thermal  flow. 

Of  course,  for  the  problems  considered  in  this  report,  neither  the  Reynolds  number  nor  the 
Rayleigh  number  can  be  determined  a-priori ;  the  convective  velocity  term  in  the  Reynolds 
number  and  the  temperature  difference  term  in  the  Rayleigh  number  must  be  determined  as 
solutions  to  the  problem. 

2.3  SIMULATION  PARAMETERS: 


The  simulation  results  for  conductive  and  convective  heat  exchange  are  obtained  using 
COMSOL  Software,  in  which  finite-element  techniques  are  used  with  the  Galerkin  method 
[Reddy  and  Gartling,  2001;  Zimmerman,  2006],  Computational  results  are  given  in  the  next  few 
chapters:  Results  for  a  two-dimensional  problem  are  given  in  chapter  3.  Results  for  three- 
dimensional  problems  are  given  in  chapter  4.  Results  for  the  three-dimensional  problem  I,  in 
which  all  boundaries  are  thermally  insulated,  are  given  in  section  4.1 .  Results  for  the  three- 
dimensional  problem  II,  in  which  heat  transfer  takes  place  at  the  front  surface  and  all  other 
surfaces  are  thermally  insulated,  are  given  in  section  4.2.  Free-surface  affects  are  discussed  in 
chapter  5.  Evaporative  effects  are  discussed  in  section  5.1,  and  Marangoni  effects,  involving 
thermally-induced  surface  tension,  are  discussed  in  section  5.2. 
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In  the  simulations  described  in  chapters  3  and  4  of  this  report,  the  following  assumptions  and 
approximations  are  made: 


1 .  Light  scattering  is  negligible  compared  to  absorption. 

2.  The  Boussinesq  approximation  is  used  with  the  incompressibility  condition  for  the  fluid. 

3.  Non-slip  boundary  conditions  (with  zero  velocity)  apply  to  all  surfaces. 

4.  Viscous  dissipative  losses  are  negligible. 

5.  Temperature  variations  of  r\,  k,  and  C  are  not  considered.  (First  Approximation)  The 
parameters  representing  the  dynamic  viscosity,  thermal  conductivity  and  specific  heat  of 
the  liquid  are  considered  as  constants  with  the  values  listed  in  table  1 . 

The  parameters  defining  all  simulations  of  chapters  3  and  4  are  given  in  table  1 . 


Table  1.  Parameter  values  used  for  all  two-dimensional  and  three-dimensional  heat  conduction/convection 
problems  [COMSOL  Lib,  2007] 


Parameter  name 

Symbol 

Value 

Units 

dynamic  viscosity 

n 

! 

o 

kg/(m.s) 

thermal  conductivity 

k 

0.6114 

W/(m.K) 

specific  heat 

C 

4200 

J/(kg.K) 

incident  laser  power 

P 

1 

W 

e_1  beam  radius 

Rw 

1.5(10“3) 

m 

absorption  coefficient 

Ma 

150 

nr1 

Density  of  water 

p  (**) 

1000 

Kg/m3 

**  According  to  the  Boussinesq  approximation,  p  is  constant,  with  the  exception  of  the  buoyant 
body-force  term  given  by  equation  (6). 


To  approximate  the  Rayleigh  number  of  water  for  small  temperature  variations,  AT,  around  295 
K,  the  coefficient  of  thermal  expansion,  a,  of  3(10-4)  [K'1],  and  g  equal  to  9.81  m/sec2  are  used. 
For  the  model  described  in  this  chapter,  the  order  of  the  distance  over  which  temperature 
variations  are  significant,  is  taken  as  dj  =  10"2  meters.  These  values  and  the  values  of  the 
parameters  given  in  table  1  are  substituted  into  equations  (8)  and  (9)  to  obtain  an  approximation 
of  the  Rayleigh  number  in  terms  of  the  temperature  variation,  AT.  The  approximate  result  is 

Ra  «  2.7(104)  AT,  (10) 

where  AT  is  given  in  units  K.  As  indicated  by  equation  (10),  large  values  of  the  Rayleigh  number 
exist  for  small  temperature  variations.  This  suggests  the  presence  of  convective  flow  for  the 
temperature  variations  within  the  size  of  the  container  considered  in  this  report. 
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3.  SIMULATION  OF  HEAT  TRANSPORT  WHEN  A  HORIZONTAL  LASER 
BEAM  IS  NORMAL  TO  A  TWO-DIMENSIONAL  PLANE  SURFACE: 

For  illustrative  purposes,  a  two-dimensional  problem  is  first  considered.  This  problem  involves 
the  simulation  of  heat  transport  when  the  heat  source  consists  of  a  horizontal  laser  beam  that  is 
normal  to  and  at  the  center  of  a  two-dimensional  plane  surface.  The  heat  source  is  represented  by 

equations  (4)  and  (5)  where  the  term  describing  exponential  decay  in  z  (the  e  /UaZ  term  in 
equation  (5)  is  set  equal  to  one.  For  this  problem,  non-slip  (zero-velocity)  fluid-dynamic 
boundary  conditions  and  thermally-insulated  thermal  boundary  conditions  exist  on  all  surfaces. 
The  initial  temperature  is  293.15  K,  and  the  initial  velocity  everywhere  is  zero. 

The  simulated  results  are  obtained  using  Comsol  finite-element  software  with  24,656  degrees  of 
freedom;  1,925  mesh  points;  and  3,728  triangular  elements. 

The  geometry  is  shown  in  Figure  2.  The  temperature  profiles  at  times  equal  to  one  and  three 
seconds  are  shown  in  figures  3  and  4  respectively.  The  horizontal  centerlines  are  drawn  in  these 
figures  to  demonstrate  the  vertical  asymmetry  of  the  temperature  profile  and  to  thereby  illustrate 
the  vertical  convective  effects  of  the  heat  transfer.  The  vertical  centerlines  are  drawn  to  illustrate 
the  horizontal  symmetry  of  the  temperature  profile.  A  convective  flow  pattern  is  shown  to  exist 
at  time  equal  to  three  seconds. 

Convective  velocity  vectors  are  shown  at  time  equals  three  seconds  in  figure  5.  Note  the  closed 
paths  formed  by  the  flow  vectors;  this  property  is  consistent  with  the  solenoidal  properties  of  the 
velocity  vector,  as  indicated  by  equation  (2). 

Temperature  profiles,  as  functions  of  y,  along  the  central  vertical  axis  (where  x=0.01m)  at 
various  times  are  shown  in  figure  6.  Note  that  the  temperature  profile  is  almost  symmetric  about 
the  central  vertical  position  (y  =  0.01  m)  at  time  equal  to  0.3  seconds.  However,  as  time 
increases,  the  vertical  asymmetry  of  the  temperature  profile  about  the  central  horizontal  position 
increases.  This  asymmetry  illustrates  the  occurrence  of  convective  heat  transfer. 

Temperature  profiles,  as  functions  of  x,  along  the  central  horizontal  axis  (where  y  =  0.01m)  at 
various  times  are  shown  in  figure  7.  These  plots  are  used  to  check  the  results  for  the  expected 
left-right  symmetry  properties  of  the  solutions. 
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y(m) 


X  (m) 


Figure  2.  The  geometry  for  the  two-dimensional  problem  is  shown.  The  laser  beam  (with  its  e-1  radius)  is  illustrated 
by  the  solid  circle.  All  boundaries  are  thermally  insulated,  and  the  fluid  velocity  at  all  boundaries  is  zero  (the  non¬ 
slip  boundary  condition  applies). 


Figure  3.  Results  for  the  two-dimensional  problem  are  shown.  The  temperature  profile  is  given  at  time  equals  one 
second.  Minimum  temperature  equals  293  K;  maximum  temperature  equals  294.6  K.  The  horizontal  centerline  is  drawn 
to  indicate  the  slight  vertical  asymmetry  of  the  temperature  distribution;  this  illustrates  an  early  stage  of  convective  flow. 
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Figure  4.  Results  for  the  two-dimensional  problem  are  shown.  The  temperature  profile  at  time  equals  three  seconds 
is  given.  Minimum  temperature  equals  293  K;  maximum  temperature  equals  299  K.  The  horizontal  centerline  is 
drawn  to  indicate  the  vertical  asymmetry  of  the  temperature  distribution;  this  illustrates  effects  of  convective  flow. 


Figure  5.  Results  for  the  two-dimensional  problem  are  given.  Convective  fluid  velocity  vectors  at  time  t  =  3s  are 
shown.  Note  the  convective  flow  profile.  The  maximum  amplitude  of  the  velocity  is  4.4  mm/sec. 
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(where  x  =  0.01  m)  are  shown  at  various  times.  The  asymmetry  (about  the  center-y-position  at  y  =  0.01  m)  of  the 
temperature  in  the  y-direction  is  indicative  of  convective  flow.  This  asymmetric  pattern  is  shown  to  increase  with 
time. 


x  [  m  ] 

Figure  7.  Results  for  the  two-dimensional  problem  are  given.  Plots  are  shown  of  temperature  profiles  along  the 
central  horizontal  axis  (where  y  =  0.01  m)  at  various  times.  These  results  are  presented  to  demonstrate  the  extent  of 
the  right-left  symmetry  (about  the  central-x  position  at  x  =  0.01  m)  of  the  result. 
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4.  SIMULATION  OF  HEAT  TRANSPORT  WHEN  A  HORIZONTAL  LASER 

BEAM  IS  NORMAL  TO  A  PLANE  SURFACE  OF  A  THREE-DIMENSIONAL 
RECTANGULAR  CONTAINER 

4.1  THREE-DIMENSIONAL  PROBLEM  1:  SIMULATION  APPLIED  TO 
THERMAL  BOUNDARY  CONDITIONS  IN  WHICH  ALL  SURFACES 
ARE  THERMALLY  INSULATED 

For  the  three-dimensional  problem  1  considered  in  this  chapter,  thermally-insulated  thermal 
boundary  conditions  are  specified  on  all  surfaces.  Also,  non-slip  fluid-dynamic  boundary 
conditions  are  specified  on  all  surfaces.  The  x-y  faces  are  perpendicular  to  the  laser  beam.  The 
laser  beam  enters  perpendicularly  to  the  front  face  at  z  =  0,  and  the  laser  penetrates  with 
exponential  decay  in  the  z-direction  (Beer’s  Law).  As  illustrated  in  figure  8,  the  faces 
perpendicular  to  the  laser  beam  are  specified  with  dimensions  0  <  x  <  0.02  m,  and  0  <  y  <  0.02 
m  and  the  penetration  distance  extends  from  z  =  0  to  z  =  0.04  m.  The  initial  temperature  is 
293. 15K,  and  everywhere  the  initial  velocity  is  zero. 


Figure  8.  The  geometry  for  the  three-dimensional  problems  1  and  2  is  illustrated.  The  arrow  shows  the  direction  of 
propagation  and  attenuation  of  the  laser  beam.  The  coordinates  are  shown  in  meters,  where 

0  <  x  <  0.02  m,  0  <  y  <  0.02  m,  and  0  <  z  <  0.04m. 

The  simulated  results  are  obtained  using  Comsol  finite-element  software  with  27,793  degrees  of 
freedom;  969  mesh  points;  and  4310  tetrahedral  elements. 

The  temperature  profile  at  the  front  face  (where  z  =  0)  at  time  equals  four  seconds  is  shown  in 
figure  9.  A  horizontal  centerline  is  drawn  to  illustrate  the  vertical  asymmetry  of  the  temperature 
profile.  The  temperature  profile  through  a  central  y-z  plane  (where  x  =  0.01m)  at  time  equals  4 
seconds  is  shown  in  figure  10. 
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Temperature  profiles,  as  functions  of  y,  along  the  vertical  central  axis  (where  x  =  0.01m)  at  the 
front  face  (z  =  0)  at  various  times  are  shown  in  figure  1 1 .  The  temperature  profile  is  almost 
vertically  symmetric  (about  the  central  y  position  at  y  =  0.01m)  at  time  equals  one  second. 
However,  the  vertical  asymmetry  of  the  temperature  profile  is  shown  to  increase  with  increasing 
time. 

Temperature  profiles,  as  functions  of  y,  along  the  central  vertical  axis  (where  x  =  0.01m)  at  time 
equals  three  seconds  are  shown  for  various  values  of  z  in  figure  12.  Note  the  asymmetry  in  the  y- 
direction.  The  vertical  asymmetry  of  the  temperature  profiles,  as  shown  in  figures  9  to  12,  are 
indicative  of  thermal  convective  effects. 

Temperature  profiles,  as  functions  of  x,  along  the  central  horizontal  axis  (where  y  =  0.01m)  at 
time  equals  three  seconds  are  shown  for  various  values  of  z  in  figure  13.  These  results  are 
plotted  as  a  check  to  determine  the  horizontal  symmetry  of  the  solution.  As  expected,  these 
results  show  that  the  solutions  are  close  to  being  symmetric  in  the  horizontal  direction  about  the 
central  horizontal  position  (at  x  =  0.01m).  The  small  deviation  from  perfect  symmetry  in  the  x- 
direction  is  believed  to  be  caused  by  a  small  computational  inaccuracy. 


\t'- 


Figure  9.  Results  for  the  three-dimensional  problem  1  are  given.  The  temperature  profile  at  the  front  face  is  shown 
at  t  =  4s.  A  horizontal  centerline  is  drawn  to  illustrative  vertical  asymmetry  caused  by  convection.  A  vertical 
centerline  is  drawn  to  evaluate  the  expected  symmetry  in  the  horizontal  direction.  Maximum  temperature  equals 
302.7  K;  minimum  temperature  equals  293.15  K 
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Figure  10.  Results  for  the  three-dimensional  problem  1  are  given.  The  temperature  profile  through  a  central  y-z 
plane  (where  x  =  0.01m)  at  t  =  4s  is  presented.  Note  the  vertical  asymmetry  of  the  temperature  profile.  Maximum 
temperature  equals  302. 8K;  minimum  temperature  equals  293.15  K. 


Temperature  [K1 
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Figure  11.  Results  for  the  three-dimensional  problem  1  are  given.  Temperature  profiles  along  the  vertical  central 
axis  (where  x  =  0.01m)  at  the  face  (z  =  0)  at  various  times  are  presented.  Note  the  increase  in  the  vertical  asymmetry 
of  the  temperature  about  the  central  y  position  (at  y  =  0.01  m)  with  increasing  time.  The  dip  in  temperature  below 
293K  is  caused  by  a  small  computational  inaccuracy. 
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Figure  12.  Results  for  the  three-dimensional  problem  1  are  given.  Temperature  profiles  as  a  function  of  y  along  the 
central  vertical  axis  (x  =  0.01m)  are  shown  for  time  equals  three  seconds.  Results  are  presented  for  various  values  of 
z.  Note  the  asymmetry  in  the  y-direction;  this  illustrates  effects  of  vertical  convection.  The  dip  in  temperature  below 
293K  is  caused  by  a  small  computational  inaccuracy. 


Temperature  [K] 


Figure  13.  Results  for  the  three-dimensional  problem  1  are  given.  Temperature  profiles  as  a  function  of  x  along  the 
central  horizontal  axis  (y  =  0.01)  at  time  equals  three  seconds  are  shown.  Results  are  presented  for  various  values  of 
z.  These  results  were  plotted  as  a  check  for  the  evaluation  of  the  horizontal  symmetry  of  the  solution. 
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4.1.1  ENERGY-CONSERVATION  CHECKS  OF  THE  SOLUTIONS  OF 
THREE-DIMENSIONAL  PROBLEM  1 

In  three-dimensional  probleml,  all  surfaces  are  specified  to  be  rigid  and  thermally  insulated; 
therefore,  energy  conservation  is  easily  accounted  for  within  the  enclosed  volume  of  the 
container.  Energy  conservation  principles,  involving  the  first  law  of  thermodynamics,  can  be 
used  as  a  means  of  checking  the  results  of  this  problem  .  Since  no  external  work  is  performed 
by  the  fixed  surfaces,  and  there  is  no  heat  loss  from  the  thermally  insulated  surfaces,  according 
to  the  first  law  of  thermodynamics, 

Energy  Input  =  Internal  Energy  Increase  (for  any  time  duration,  td), 

When  the  laser  power  input  is  constant,  during  a  laser  time  duration,  td ,  the  energy  input  is  given 
as 

Energy  Input  =  td  JJJ  Qs  (x,  y,  z)  dV ,  (11) 

Vol 

integrated  over  the  volume  of  the  container. 

The  Energy  Input  [J]  can  be  determined  by  substituting  equations  (4)  and  (5),  with  the 
parameters  given  in  table  1,  into  the  integral  expression  (11).  Evaluating  the  integral  over  the 
rectangular  volume  of  the  container,  we  get  the  result  for  a  laser  time  duration  td  and  input  laser 
power,  P : 


Ei  =  Ptd  erf 


f 


L 


X 


2  R 


erf 


W  ) 


f  Ly 
\2RW  j 


_e~PaLZ 


(12) 


where  erf )  denotes  the  error  function.  Substituting  fia  =  150  m"1,  P  equals  1W,  Rw  equals  1.5 
(10"  )  m  with  the  container  dimensions:  Lx  equals  0.02m,  LY  equals  0.02m,  and  Lz  equals  0.04 
m,  we  get  the  following  result  the  energy  input  during  time  duration,  td  : 

E{  =  0.9975  td  (12a) 


The  internal  energy  increase  [J]  of  the  liquid  during  the  laser  time  duration  td  [s]  can  be 
calculated  as 


Internal  Energy  Increase  =  pC  JJJ  [T (x,  y,  z)  -  7q  ]dV 

Vol 


(13) 


For  the  small  convective  velocities  involved,  viscous  dissipation  of  energy  can  be  neglected. 
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where  T(x,y,z)  is  the  calculated  temperature  at  time  tj,  and  To  is  the  specified  initial 
temperature  (To  =  293. 15K).  The  integration  is  over  the  volume  of  the  container. 

The  internal  energy  increase  is  determined  by  substituting  p  =  1000  kg/m  ,  C  equal  to  4200 
J/(kg»K)  and  the  initial  condition,  T0  equal  to  293. 15K,  into  the  integral  expression  (13)  .  Here 
the  integral  is  evaluated  using  the  post-processing  capability  of  COMSOL  software,  where  a 
numerical  integration  is  performed  using  the  evaluated  results  for  T(x,y,z). 

Calculated  results  for  input  energy,  internal  energy  increase,  and  the  percent  difference  between 
the  energy  input,  Eit  and  the  internal  energy  increase  in  are  shown  in  terms  of  the  specified  time 
durations  in  table  2.  These  results  appear  to  show  reasonable  agreements  between  the  calculated 
energy  input  and  the  calculated  internal  energy  increase,  especially  for  laser  time  durations 
greater  than  one  second. 


Table  2:  Results  of  the  energy  conservation  check.  Laser  time  duration,  input  energy,  internal  energy  increase,  and 
percent  difference  between  the  input  energy  and  the  internal  energy  increase  of  the  water  in  the  closed  container. 


td  (sec) 

Input  Energy 
(Joules) 

Internal  Energy  Increase 
(Joules) 

Percent  Difference 

1 

0.998 

0.896 

10.2 

2 

1.996 

1.854 

7.1 

3 

2.794 

2.919 

4.6 

4 

3.992 

3.897 

3.4 

5 

4.990 

4.851 

2.9 

6 

5.988 

5.771 

3.6 

7 

6.986 

6.690 

4.2 

8 

7.984 

7.610 

4.7 

4.2  THREE-DIMENSIONAL  PROBLEM  2:  SIMULATION  APPLIED  TO 
THERMAL  BOUNDARY  CONDITIONS  WHERE  HEAT  TRANSFER 
OCCURS  AT  THE  FRONT  SURFACE  (AT  Z  =  0)  AND  ALL  OTHER 
SURFACES  ARE  THERMALLY  INSULATED 

The  problem  considered  in  this  chapter  is  identical  to  the  three-dimensional  problem  1,  treated  in 
section  4.1,  with  the  exception  that  here  heat  transfer  between  the  container  and  its  surroundings 
occurs  at  the  front  surface  where  z  =  0.  Thermally-insulated  thermal  boundary  conditions  are 
specified  on  the  other  surfaces.  Also,  similar  to  problem  1,  non-slip  fluid-dynamic  boundary 
conditions  are  specified  on  all  surfaces.  The  x-y  faces  are  perpendicular  to  the  laser  beam.  The 
laser  beam  enters  perpendicularly  at  the  center  of  the  front  face  at  z  =  0,  and  the  laser  penetrates 
with  exponential  decay  in  the  z-direction.  As  illustrated  in  figure  8,  the  faces  perpendicular  to  the 
laser  beam  have  specified  dimensions  0  <  x  <  0.02  m,  and  0<y<  0  .02  m  and  the  penetration 
distance  extends  from  z  =  0  to  z  =  0.04  m.  The  parameters  specified  in  table  1  apply. 
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The  heat  transfer  at  the  z=0  surface  is  described  by  Newton’s  Law  of  Cooling  [Incropera,  et  al, 
2007;  Wolfram  Research,  2006].  The  initial  temperature  in  the  container  is  293. 15K,  and  the 
initial  velocity  is  zero  everywhere.  Convective  heat  exchange  occurs  between  the  surface  at  z=0 
and  the  surroundings,  which  are  at  a  fixed  temperature  of  293. 15K.  The  thermal  boundary 
condition  at  the  z=0  surface  is  specified  with  a  heat-transfer  rate  of  QSurf  [W/m  ]: 

Qsurf  = |z=0- 293.15),  (14) 

where  T|z=o  is  the  temperature  of  the  fluid  at  z  =  0+,  and  the  parameter  he  is  specified  as  he  = 
15  [W/(m2»K]  **  The  boundary  value,  T|z=o  ,  interactively  computed  at  z  =  0,  is  dependent 
upon  the  coordinates  x  and  y  and  time  t.  The  sign  convention  is  that  heat  flow  into  the  fluid  at 
the  z  =  0  boundary  is  positive;  this  explains  the  minus  sign  in  equation  (14). 

The  simulated  results  are  obtained  using  Comsol  finite-element  software  with  27,793  degrees  of 
freedom;  969  mesh  points;  and  4310  tetrahedral  elements. 

The  temperature  profile  at  the  front  face  (where  z=0),  at  time  equals  four  seconds,  is  shown  in 
figure  14.  A  horizontal  centerline  is  drawn  to  illustrative  vertical  asymmetry  (about  the  central 
horizontal  axis  at  y  =  0.01m)  due  to  vertical  convection.  A  vertical  centerline  is  drawn  to 
evaluate  the  expected  symmetry  in  the  horizontal  direction. 

The  temperature  profile  through  the  central  vertical  y-z  plane  (where  x  =  0.01m)  at  time  equals 
four  seconds  is  shown  in  figure  15.  Here,  the  horizontal  line  is  drawn  to  show  vertical 
asymmetry  in  temperature  due  to  thermal  convection. 

Figure  16  shows  the  calculated  temperature  profiles  along  the  vertical  central  axis  (at  x=0.01m) 
at  the  front  surface  (z  =  0)  at  various  times.  Figure  17  shows  the  calculated  temperature  profiles, 
at  various  values  of  z,  as  a  function  of  y  along  the  central  vertical  axis  (where  x=0.01)  at  time 
equals  four  seconds.  The  vertical  center,  at  y=0.01m,  is  drawn  on  figures  16  and  17  to  illustrate 
the  vertical  asymmetry  of  the  calculated  temperature  profiles. 

Temperature  profiles  as  a  function  of  x  along  the  central  horizontal  axis  (where  y  =  0.01m)  at 
time  equals  four  seconds  are  shown  for  various  values  of  z  in  figure  18.  These  results  are  plotted 
as  a  check  to  evaluate  the  horizontal  symmetry  of  the  solution.  The  results  demonstrate  the  near 
horizontal  symmetry  (about  the  central  position  at  x  =  0.01  m)  of  the  calculated  temperature 
profiles.  Slight  deviations  from  perfect  horizontal  symmetry  can  be  attributed  to  small 
inaccuracies  in  the  computational  method. 


The  value  specified  by  the  parameter  he  is  arbitrary.  The  value  of  he=15  corresponds  to  the  heat-transfer 
coefficient  from  skin  to  air.  The  value  of  he  for  the  air/container  interface  is  dependent  upon  the  thermal-insulation 
properties  of  the  wall  of  the  container  and  its  surroundings.  Since  this  value  is  unknown,  a  value  of  he=15  is  chosen 
here. 
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Figure  14.  Results  for  three-dimensional  problem  2  are  given.  The  temperature  profile  at  the  front  face  at  time 
equals  four  seconds  is  shown.  A  horizontal  centerline  is  drawn  to  illustrate  vertical  asymmetry  (about  the  central 
horizontal  axis)  due  to  thermal  convection.  A  vertical  centerline  is  drawn  for  the  evaluation  of  the  expected 
symmetry  in  the  horizontal  direction.  Maximum  temperature  equals  296.2  K;  minimum  temperature  equals 
293. 15K 


Mta.iSA.lS9 


Figure  15.  Results  for  three-dimensional  problem  2  are  given.  The  temperature  profile  through  the  central  vertical 
plane  (at  x  =  0.01m)  at  time  equals  four  seconds  is  presented.  The  horizontal  line  is  drawn  to  show  vertical 
asymmetry  due  to  convection.  Maximum  temperature  equals  296. 2K;  Minimum  temperature  equals  293. 15K. 
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Figure  16.  Results  for  three-dimensional  problem  2  are  given.  The  temperature  profiles  along  the  vertical  central 
axis  (at  x  =  0.01m)  at  the  front  face  (where  z  =  0)  at  various  times  are  presented 


Figure  17.  Resuus  ror  mree-aimensionai  prooiem  z  are  given,  me  temperature  promes,  as  a  iunction  or  y,  along 
the  central  vertical  axis  (where  x  =  0.01)  at  time  equals  four  seconds  are  presented.  Results  are  shown  for  various 
values  of  z.  Note  the  asymmetry  in  the  y-direction;  this  illustrates  effects  of  vertical  convection. 
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Temperature  [K] 


Figure  18.  Results  for  three-dimensional  problem  2  are  given.  Temperature  profiles  as  a  function  of  x  along  the 
central  horizontal  axis  (y  =  0.01)  at  time  equals  four  seconds  are  presented.  Results  are  shown  for  various  values 
of  z.  These  results  are  plotted  as  a  check  for  the  evaluation  of  the  horizontal  symmetry  of  the  solution. 
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5. 


FREE-SURFACE  EFFECTS 


The  problems  described  in  chapter  4  involve  heat  transfer  in  closed  containers.  For  the  problem 
considered  in  section  4.1,  all  boundaries  are  thermally  insulated.  For  the  problem  considered  in 
section  4.2,  convective  heat  transfer  occurs  between  the  front  boundary  (at  z  =  0)  and  its 
surroundings.  The  fluid-dynamic  boundary  conditions  for  problems  described  in  sections  4.1  and 
4.2  involve  zero  velocity  at  the  walls  (the  non-slip  conditions  for  fixed  surfaces). 

In  this  chapter,  the  effects  of  free-surface  boundary  conditions  are  discussed.  These  surface 
effects  include  convective  and  evaporative  heat  transfer  and  surface-tension  induced  by 
temperature  gradients  on  the  surface  (the  Marangoni  Effect). 

5.1  EVAPORATIVE  AND  CONVECTIVE  HEAT  TRANSFER  AT  A  FREE 
SURFACE: 

Two  mechanisms  of  energy  exchange  can  take  place  at  the  air-liquid  surface,  at  z  =  0.  One 
mechanism  involves  the  convective  heat  transfer  between  the  liquid  at  density  Pl  (kg/m3)  and  the 
vapor  at  density  pv  (kg/m  ).  The  other  mechanism  involves  heat  loss  from  water  evaporation  at  z 
=  0  [Tilbas  and  Sami,  1998],  Thus,  the  thermal  boundary  condition  at  z  =  0  (at  the  air-water 
interface)  is 


k8T 

dZ 


Z=0 


+  Qsurf  • 


(15) 


The  heat  evaporation  term,  (2vap 


[W/m2],  at  z  =  0  is 


*m 


(/^v(sat) 


(16) 


Here  hfg  is  the  parameter  representing  enthalpy  change  due  to  evaporation  (J/kg),  hm  is  the  mass- 
transport  coefficient  (m/s),  pv(sat)  is  the  mass  density  of  the  water  vapor  at  the  water  surface,  and 
Pv(inf)  is  the  mass  density  of  the  water  vapor  at  the  room  temperature.  The  convective  heat 

transfer  term  at  z=0,  Qsurf ,  is  given  by  equation  (14). 


Consider  a  cylindrical  column  of  water  with  a  laser  beam  heat  source  and  with  heat  transfer  due 
to  evaporation  and  convection  at  the  top  (z  =  0)  surface  as  shown  in  figure  19. 
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Cylindrically  Symmetric 
Laser  Beam 


Zero  heat  flux 

Figure  19.  The  geometry  for  the  simulation  involving  the  evaporative  and  convective  heat  transfer  at  the  z  =  0 
surface  is  illustrated.  The  height  of  the  cylinder  equals  1  cm,  and  the  radius  of  the  cylinder  equals  0.5  cm.  The  laser 
beam  enters  perpendicularly  at  the  top  (at  z  =  0)  surface.  The  problem  is  circularly  symmetric. 


The  heat  conduction  equation  is 

pc--V*(kVT)=Qs.  (i7) 

For  simplicity,  heat  convection  within  the  fluid  is  not  considered  in  this  chapter. 

The  heat  source  term,  Qg,  is  described  by  equations  (4)  and  (5)  (with  appropriate  alterations 

corresponding  to  the  cylindrical  geometry  for  the  problem  solved  in  this  chapter).  The  problem 
is  solved  in  cylindrical  coordinates  with  axial  symmetry. 

The  problem  is  solved  using  COMSOL  software.  As  evaporation  takes  place,  the  water  layer  is 
depleted  at  the  evaporative  surface  at  z=0.  To  account  for  the  depleted  water  layer,  the  problem 
is  solved  with  moving  coordinates  (moving  boundaries).  That  is,  at  z=0,  the  mass  depletion  per 
unit  surface  area  per  unit  time  (kg/(m2s) )  is  calculated  as 


M  ~  hm  (Pv(sat)  Pw(inf)  )• 


(18) 
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The  surface  velocity  of  the  moving  boundary  at  z  =  0  is 


VZ|z=0 


M 

pl’ 


where  pL  is  the  density  of  the  liquid  water  at  the  evaporating  surface. 


(19) 


The  parameters  describing  this  problem  are  shown  in  Table  3a,  and  the  relations  used  are  shown 
in  Table  3b.  The  simulated  results  are  obtained  using  Comsol  finite-element  software  with  3,271 
degrees  of  freedom;  277  mesh  points;  and  406  triangular  elements. 


Results  showing  the  temperature  profiles  at  time  equals  eight  seconds  are  shown  in  figure  20. 
Alterations  in  the  shape  of  the  top  surface  were  small  and  were  not  apparent  in  figure  20  for  the 
incident  laser  power  (1W)  considered.  Alterations  in  the  shape  of  the  top  surface  would  be 
observed  if  higher  laser  power  were  used  in  the  simulation. 


Table  3a.  Parameters  describing  the  free-surface  problem 


Parameters  Value  Description 


Tinf 

293.15 

Room  Temp  (K) 

he 

15 

Heat  convection  coeff  (W/mA2K) 

hfq 

2 . 35e6 

phase-chanqe  enthalpy  (J/kqK) 

Rh 

0.50 

Relative  humidity 

rho  v  sat  inf 

0.0173 

mass  density  of  saturated  water 
vapor  at  room  temperature 

P  beam 

1 

Laser  power  (w) 

mu 

150 

absorption  coeff  (mA-l) 

Rw  (R  waist) 

0.0015 

beam  waist  size  (m) 

Rho  L 

1,000 

density  of  liquid  water  (kg/mA3) 

T  room 

293.15 

room  temp  (K) 
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Table  3b.  Relations  between  the  expressions  describing  the  evaporation  problem 


Expressions 


Term 

Expression 

Description 

hm 

he*  (le-6*Tcc+0 .000  9) 

Convection  mass 

transfer  coeff 
(m/s) 

Qsource 

2.*P  beam/ (pi*r  waist A2 ) *exp  (- 
2 . *r A2/  (r  waist A2 )  ) *exp (-mu*z ) 

Laser  source  power 
density  (W/mA3) 

Tcc 

T-273 

Conversion  from 
cent.  To  Kelvin 

Temp 

pho  v  sat  surf 

le-3* (4e-6*TccA4-6e- 
5*TccA3+0 . 0196e- 
3*TccA2+0 . 1534 *Tcc+6 . 1098) 

Density  of 
saturated  water 
vapor  at  z=0 
surface 

Qvap 

hfg*hm* (pho  v  sat  surf- 
rho  v  sat  inf) 

Evap  surface  heat 
losses  (W/mA2) 

Qconv  surf 

he* (T-Tinf ) 

Conv.  surface 
losses  (W/mA2) 
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Incident  laser  beam 


:< —  0.5  cm  — ► 

i 

Center 

Figure  20.  The  simulated  temperature  profile  in  a  cylindrical  container  of  water  is  shown.  Evaporative  cooling 
and  convective  heart  exchange  take  place  at  the  top  boundary,  where  z  equals  zero.  Results  for  time  equals  eight 
seconds  are  presented.  Laser  input  power  equals  1 W.  Maximum  temperature  equals  3 12K;  Minimum  temperature 
equals  293. 15K 


5.2  SURFACE  TENSION  AT  A  FREE  SURFACE:  THE  MARANGONI 
EFFECT 

Surface  tension  at  the  surface  of  a  liquid  may  be  caused  by  radius-of-curvature  variations  along 
the  surface  and/or  by  temperature  gradients  along  the  surface  [Pimputkar  and  Ostrach.,  1980], 

The  production  of  surface  tension  by  variations  of  the  radius  of  curvature  along  a  surface  is 
described  by  Laplace’s  Law  [Landau  and  Lifshitz,  1959].  This  effect  is  commonly  small  for  the 
temperature  range  considered  in  this  report. 

The  production  of  surface  tension  by  temperature  gradients  along  a  free  surface  is  described  as 
the  Marangoni  Effect  [Levich,  1962;  Levich  and  Krylow,  1969;  Longtin  et  al,  1999].  The  surface 
boundary  condition  at  the  liquid-air  interface  is  given  in  terms  of  the  surface  stresses  that  are 
produced  by  Marangoni  convection;  these  stresses  at  the  liquid-air  interface  are  given  as 

du  r  dT 

r] — —  =  y - 
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(20) 


where  x  is  the  direction  tangent  to  the  surface  and  y  is  the  direction  normal  to  (and  into)  the 
surface.  The  parameter  y  is  equal  to  the  derivative  of  the  surface  tension  with  respect  to 
temperature.  A  plot  of  the  surface  tension  vs.  temperature  of  water  is  shown  in  figure  21.  The 
parameter  y  is  determined  by  taking  the  derivative  with  respect  to  temperature  of  the  curve 
shown  in  Figure  21. 

The  surface  tension  would  drive  fluid  flow  at  and  near  the  free  surface  [Yih,  C.S,  1968,  1969].  A 
simulation  showing  the  flow  induced  by  Marangoni  convection  is  given  in  the  COMSOL 
Multiphysics  Modeling  Guide  [2005]. 


Figure  21.  Surface  tension  vs.  temperature  of  water.  The  value  of  the  parameter  y  can  be  found  by  taking  the 
derivative  of  this  curve. 
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6. 


SUMMARY 


Mathematical  models  and  their  simulated  results  are  presented  to  describe  heat  transfer  from  a 
penetrating  laser  beam  in  a  rectangular  container  of  water.  Beer’s  law  is  used  as  an 
approximation  of  the  absorption  loss  of  the  laser  beam,  and  scattering  effects  are  neglected.  Both 
convective  and  conductive  heat  transfer  are  considered.  COMSOL  software  is  used  to  find  finite- 
element  simulation  results  for  convective  flow  velocity  and  temperature  changes  within  the 
container.  Both  simulated  results  and  dimensional  analysis  indicate  that  convective  flow  is 
important. 

The  following  suggestions  are  made  concerning  the  manner  in  which  this  work  can  be  extended 
in  future  investigations: 

1 .  The  simulated  temperature  changes  could  be  coupled  with  calculated  optical  effects  in 
which  changes  of  refractive  index  with  temperature  can  be  taken  into  account.  The 
importance  of  thermal  convection  on  the  optical  measurements  has  been  suggested  in 
experiments  by  Vincelette  et  al  [2007], 

2.  The  simulated  temperature  changes  should  be  verified  experimentally.  For  this  purpose, 
the  temperature  mapping  methods  of  Maswadi  et  al  [2004]  could  be  used. 

3.  The  work  can  be  extended  to  incorporate  surface  phenomena,  as  described  in  chapter  5, 
with  the  heat  conduction/convection  simulations  presented  in  chapter  4. This  work  can  be 
further  extended  to  include  the  simulation  of  temperature  changes  and  the  associated 
changes  in  the  optical  properties  of  a  tear  film  at  the  surface  of  the  cornea  of  the  eye. 
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